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A computational fluid dynamics methodology is used to study the salient flow features around 
the breasts of a human figure and to describe the aerodynamic differences imparted by their 
geometric presence. Two models are proposed for examination: a 3-dimensional reference based 
on a character design with a significantly buxom figure and a modification of this design where 
the breast size is reduced significantly. The two models are tested at speeds ranging from 1 to 30 
m-s A -l using Reynolds-averaged Navier Stokes (RANS). Drag, lift, and skin friction forces, 
along with turbulence kinetic energy (TKE), are investigated and compared between the differ¬ 
ent models. The present results are expected to provide useful information on the validity of the 
statement, “Flat is Justice” in terms of an aerodynamic standpoint. In addition to this, the results 
can offer worthwhile data investigating the anthropometrical presence of large breasts on sport 
aerodynamics. 


1. Introduction 

The aerodynamics of the human form has been an area of valuable 
research in various aspects of sports and competition. Air resistance 
(hereinafter referred to as “drag”) is a concerning factor in many time- 
based trials, and enhancing potential efficiency can be done through the 
elucidation of the flow around the human figure. Studies concerning the 
drag of the human body using wind tunnels can be found dating back to 
the 1920s [1]. A small sampling of subsequent studies exploring the 
effect of drag covers areas such as running [2], cycling [3], skiing [4], 
and skating [5], all of which reinforces the relevance of aerodynamic 
investigation on the human shape in regards to performance. 

In many of such studies, the authors seek to investigate the effect of 
positioning in relation to drag [6], and some utilize numerous subjects 
of differing anthropometric proportions to describe a generalized result 
on such positioning [7, 8]. Hitherto, none within the author’s investiga¬ 
tions has described the effect of specific physiological features on aero¬ 
dynamic performance in great detail. Stemming from certain internet 
communities and pertinent to the current era comes the succinct state¬ 
ment, “Flat is Justice”, which consequentially begets interesting debate 
that can reverberate and diffuse throughout media. Essentially, the 
statement describes the appreciation of flat-chested women [9], which 
posits a peculiar aspect that has yet to be fully explored in human aero¬ 
dynamics; namely, the effect of breasts in regards to drag and overall 
aerodynamic performance. 

This work is intended to contribute to the understanding of how 
large breasts can affect the dynamics of the human wake through the 
use of computational fluid dynamics (CFD) simulation tools. This pre¬ 
liminary work focuses solely on comparing the relevant effects of large 
breasts of a selected human design to that of the same design but with, 
euphemistically, “lesser tracts of land”. The following sections will 
present an overall understanding on the human wake in relation to sim¬ 
plified geometry along with engineering applications, introduce the 
chosen human geometry and models, relevant boundary conditions, the 
governing equations, and the numerical methods used to solve the equa¬ 
tions. An in-depth review on the computational uncertainty is described, 
following with extensive results and discussion, conclusions, and rec¬ 
ommendations for future work. 

1.1 Background on the Human Wake 

The human body can best be described as a bluff-body in respect to 
the flow around it. Literature on the behavior of the wakes behind bluff 
bodies indicates that the flow will be unsteady due to the turbulent tran¬ 
sition and separation of the boundary layer [10]. A simplification analo¬ 


gous to the human shape can be represented by a grouping of uniform 
circular cylinders [11] and therefore existing studies on this type of 
geometry can provide general insight into the wake region. Sumner et 
al. [12] described the wake and development of vortex structures of 
cylinders with aspect ratios (i.e. height to diameter) of 3, 5, and 9, and 
determined that a transition in vortex shedding occurs at hid = 3. An 
investigation by Okamoto and Sunabashiri [13] also supports this find¬ 
ing, adding that cylinders with an aspect ratio of 3 experience a recircu¬ 
lation region that extends four diameters downstream. Assuming the 
human form takes on a roughly large cylindrical shape near this aspect 
ratio, it is to be expected that the recirculation region will behave simi¬ 
larly and extend approximately four body widths downstream. 

A readily apparent deviation in geometry compared to studies done 
on singular cylinders is the presence of the gap between the legs. An 
extensive and comprehensive review done by Zhou and Alam [14] on 
the various arrangements of two cylinders indicate the wake structure 
falls into a multitude of regimes. In a side-by-side configuration, being 
similar to the two legs of a human, it is deduced that there are three 
primary regimes where the wake experiences proximity interference. 
When closely spaced together, the first regime shows that the cylinders 
act similarly to that of a single bluff body with a width corresponding to 
the two cylinders. When the gap width is larger than 20% of the diame¬ 
ter, each cylinder has individual wakes that strongly affect one another 
and is associated with the second regime. At gap widths exceeding ap¬ 
proximately 100-120% of the diameter, each cylinder acts as an inde¬ 
pendent body with the vortex streets being loosely influenced by one 
another. Seeing that human legs are not strictly cylinders with a fixed 
diameter but more akin to inverted tapered cylinders, the wakes behind 
the legs will likely behave in a similar fashion observed in both the first 
and second regime. With the ankles and calves being narrower and hav¬ 
ing a larger gap between them, the second regime is applicable. A tran¬ 
sition into the first regime can be expected associated with the bulkiness 
of the thighs and reduction in gap width. 

Engineering literature can also provide additional details on the 
flow characteristics around the body. Many of such studies are motivat¬ 
ed by exposure control and contaminant transport [15, 16], thermal 
issues [17, 18], and comfort prediction [19], rather than overall drag 
effects. Inherently, many of the tested flow characteristics are evaluated 
in a quiescent environment or at air velocities that are of a lower order 
compared to those found in sport-related studies. Nonetheless, these 
studies provide useful insight on the natural turbulence caused by the 
human form and the expected anatomical location of flow separation. 
Inthavong et al. [20] utilized a high speed camera to record the wake 
generation of a l/5th scaled realistic human manikin that was accelerat¬ 
ed to a velocity of ~1 m-s' 1 . From their results, it was found that the 


1 


Copyright © 2018 N. Rabino 








(a) Original reference design. (b) The modified design proposed 

Courtesy: Kyoto Animation. for comparison. 

Figure 1. Comparison of different designs for Lucoa. 



Figure 2. A perspective measurement of Lucoa in reference to a door frame 
using the Vanishing Point Tool in Adobe Photoshop. 


shoulder undergoes flow separation and produces vortices in a regular 
pattern. The hands produce a well-defined yet unstable vortex sheet that 
curls towards the centerline of the body. The head acts similarly to clas¬ 
sical sphere/cylinder cases with the addition of a trailing wake forming 
from behind the neck. The neck was found to remove the expected 
counter free shear layer that is present in cylinder studies and thus elim¬ 
inates the formation of an oscillating vortex sheet. In all, it can be said 
that the observed human wake is a highly complex and richly diverse 
system that is easily influenced by the inherent geometry used; it is 
expected that from this study, an overall summary can be presented on 
how and to what degree the previously described flow structures are 
affected by the presence of large breasts. 

2. Methodology 

2.1 Design Proposal and Model Scaling 

The use of realistic human models affords greater realization of the 
pertinent flow characteristics as they are considerably different than 
those of generalized models. Yan et al. [21] concluded that an excessive 
degree of simplification in using a manikin can affect the ability to 
achieve accurate results, and thus precludes the use of a simplified 
model for this study. However, the acquisition of a 3-D scanned human 
model with a significant bust indubitably proved difficult. The use of a 
highly unconventional approach was used to ameliorate this issue. 

The animated adaptation of Miss Kobayashi s Dragon Maid , being 
a recently popular show [22] and spawning a sizable subculture on the 
internet [23], proved suitable in terms of providing potential models. 
The dragon characters (themselves being based off of mythically and 
culturally prominent dragons) assume a human form to interact with 
other humans in this well-received [24] slice-of-life urban fantasy. A 
majority of the human forms of the female dragon characters possessed 
significant busts. However, Quetzalcoatl (referred to canonically as 
“Lucoa” and will be named as such throughout the rest of this paper) 
substantiated herself as the adaptation’s gag character by her significant 
size [25, 26], thus making her the perfect candidate in providing a suita¬ 
ble model. Being clearly the largest amongst her fellow dragons as es¬ 
tablished in Figure 1, Lucoa provides the best contrast between a large 
bust and having none at all. To provide the most direct comparison in 
regards to the effect of large breasts on the wake, a dramatic reduction 
in bust size as reflected in Figure 1(b) was proposed for use in this 
study. 


In order to obtain accurate results from the setups described later in 
this paper, it is important to have the subject in question reflect real 
world scales properly. Lucoa’s height while in human form is not given 
explicitly in any related media within the scope of the author’s research. 
Thus, Lucoa’s height must be estimated in relation to objects of which a 
reasonable measurement can readily be found. Conveniently, there is a 
scene found within Episode 6, Season 1 of the animated adaptation 
wherein Lucoa steps through a doorframe. Assuming the door is of a 
typical size 1 used for external entrances, in addition to Lucoa being 
scaled properly in the scene, we can estimate her height using a vanish¬ 
ing point technique. 

Using the door as depicted in Figure 2 to judge Lucoa’s height, it 
was determined that she stands approximately 177 cm measured to the 
top of her hat, with her horns boosting her overall figure to a height of 
182 cm. These numbers can be considered reasonable based on canoni¬ 
cal descriptions of Lucoa’s towering stature [27] compared to the aver¬ 
age height of 158 cm for a Japanese woman [28]. 

2.2 3-Dimensional Models and Geometry Analysis 

Since Lucoa is a fictional character that is commonly portrayed in 
a 2-dimensional 2 world, determining her form drag between the two 
proposed designs as described in Figure 1 requires that we add another 
dimension to her model. Conveniently enough, an available 3-D model 
of Lucoa [29] was used that would make the simulation possible. This 
MikuMikuDance 3-D model (henceforth referred to as the “Normal” 
model) was then imported into the 3-D modeling program Blender, 
scaled to the determined height as described in the previous section, 
then exported into an STL file. This STL file was then repaired using 
the built-in repair feature present in Microsoft 3D Builder due to the 
unclean geometry inherent with the model. To achieve the modified 
design (henceforth referred to as the “Flat” model), the original Miku¬ 
MikuDance model was modified using the built-in tools in Blender to 
dramatically reduce Lucoa’s breast size. The export and repair process 
remained the same as for the original model. 

As shown in Figure 3, all positions between the two models remain 
the same and left unperturbed to leave the reduction in breast size as the 
sole geometric difference to be investigated. Although the typical or¬ 
thostatic (standing) orientation of a human has the upper limbs in a 


1 A typical metric external door’s size is 926 mm wide by 2040 mm tall. 

2 Referring to the media she is portrayed in, such as printed materials and televi¬ 
sion. 
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(a) Reference (Normal) model. (b) Modified (Flat) model. 

Figure 3. 3-D representations of Lucoa to be used in CFD simulations, detailing 
(clockwise) top, side, and front views. 


relaxed position [30], the arms are left posed at a 45° adduction angle 
from the torso, as this is the default ‘A’ pose when importing the model. 
This arm position also has an advantage in this study as it potentially 
enables a more thorough analysis on the effect of breasts on the wake 
region, whereas a neutral standing posture would have the arms inter¬ 
fere with the downstream effect of the breasts. The hair is left modeled 
as solid to reduce simulation complexity and setup. While humans 
naturally lean forward against the direction of the wind to maintain 
equilibrium [31], this factor is not considered in this study as this lean¬ 
ing would change the frontal area exposed to the fluid flow and thus 
complicate comparisons against static reference models. 

Dimensionally, the bounds of the two models are similar, with the 
height and arm span being 1.82 and 1.387 meters respectively. The 
Normal model has a depth of 0.525 meters whereas the Flat model is 
only 0.414 meters. The frontal projected area, A, of both models is 
0.584 m 2 . The volumetric difference between the two is 9.19 L, indicat¬ 
ing that each breast on the Normal model has an enormous volume of 
approximately 4.6 L. The under-bust circumference of the Normal mod¬ 
el is approximately 64 cm and the bust measures 115 cm. The Flat mod¬ 
el has the same under-bust measurement whereas the bust measures 68 
cm. Attempting to match the dimensions and bust volume of the Normal 
model to existing cup sizing scales is difficult as these measurements 
are exceptionally large and exceed volumes measured in other studies 
[32]. Using the JIS L 4006:1998 [33] scale and extrapolating 3 4 cup siz¬ 
ing from the largest listed size (I-cup), the Normal model can be de¬ 
scribed as being 10 cups larger; an estimated “S65”. The Flat model is a 
stark contrast to this, where it matches a petite “AA65” size. 

The dramatic difference in bust size between the models serves to 
provide the most significant change in outcomes; it is assumed that due 
to the absurd bust size, any size smaller than the Normal model would 
have an outcome that would fall in a range between both models. 

2.3 Evaluated Metrics and Implementation 

Four metrics under investigation for this study include drag and lift 
forces (including their associated coefficients), skin friction coefficient, 
and finally, turbulence kinetic energy. To evaluate the drag coefficient, 
C D , and drag force, fy>, the following equations are used, 


Cd _ 
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where p is the fluid density, U ^ is the free stream velocity, A is the 
frontal projected area, and P is the pressure at the surface dA. r w is the 
local wall shear stress being defined as, 


with p as the dynamic viscosity, u the flow velocity along the boundary, 
and y being the height above the boundary. The value of C D is not con¬ 
stant and is dependent on Reynolds number, which is defined as, 
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where L is an arbitrary characteristic length. In this study, L is equal to 
the height of the models. 

The lift coefficient is comparable to the drag coefficient, being that 
the force is evaluated in a direction that is perpendicular to the mean 
flow direction, e.g. vertically upwards. Thus, 
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Instead of the frontal projected area, A, a reference surface area, S, is 
used. For consistent comparison however, A and S are left defined as 
being equivalent, thus A = S. This result does not affect the calculated 
forces but rather only the coefficient, and as such, the lift coefficient is 
dependent on the frontal area. 

The skin friction coefficient, C f , is evaluated in a similar manner to 
the drag coefficient since the force attributed to skin friction is a com¬ 
ponent of the profile drag, F D . Therefore, 


C f 


2 
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Analyzing the skin friction coefficient allows insight into areas where 
the boundary layer thickness changes; as turbulent flow increases, the 
thickness of the boundary layer increases, and consequently areas where 
C f transitions to larger values or experiences spikes are indicative of 
where flow separation is prevalent [34, 35]. 

Turbulence kinetic energy (TKE) signifies of the loss of kinetic en¬ 
ergy from the mean flow and represents the energy present with eddies 
in turbulent flow; it is a direct measure of the intensity of turbulence. In 
a general form quantifying the mean of turbulence normal stresses, TKE 
is defined as, 


fc = i (( U ') 2 + ( y ') 2 + ( w ') 2 ) (8) 

The exact value of TKE is calculated based on the closure of the Reyn¬ 
olds-averaged Navier-Stokes equations, which is further discussed in 
Section 3.3. 

The numerical simulations in this present work, along with the au¬ 
tomatic evaluation of the equations described in this section, were car¬ 
ried out using ANSYS FLUENT R17. The 3-D models defined in Sec¬ 
tion 2.2 were imported into FLUENT and followed the methodology as 
described in the following section. 


Ib = J (,—P cos 0 + t w sin 0) dA (2) 


3 Hair physics is beyond the scope of the author, and thus this study, due to the 
inordinate amount of computing resources and time needed to setup and simulate 
hair strands in a physically accurate fashion. 

4 In [33], each cup size is binned with every 2.5 cm deviation from the under¬ 
bust measurement starting from 7.5 cm. 


3. Computational Fluid Dynamics (CFD) Setup 
and Analysis 

3.1 Boundary Conditions 

The use of boundary conditions based on real-world environments 
enhances the overall applicability of the results stemming from the sim¬ 
ulations. It was therefore important to determine the most appropriate 
and accurate environment in which to simulate the models with. It was 
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Figure 4. Boundary conditions of the computational domain, with the inlet being 
represented in blue, outlet in red, walls in white, and symmetry in yellow. 


found that the overall location used in the animated adaptation of Miss 
Kobayashis Dragon Maid was based on the city of Koshigaya [36], 
situated in the Saitama Prefecture of Japan. A logical time of year to 
assume a person being outside without excess clothing would be some¬ 
time in the summer. Using the month of August, it was found that 
weather conditions in Koshigaya and nearby surrounding regions fea¬ 
ture averages [37] of 22.6°C for temperature, 73% for relative humidity, 
and 1005.9 hPa for local atmospheric pressure. Thus, the air density was 
calculated to be p = 1.1581 kg • m -3 and the dynamic viscosity to be 
g = 1.86847 • 10 -5 kg • m -1 • s" 1 . 

Since the human body can vary based upon the clothing worn, sur¬ 
face roughness and the effects of fabrics are parameters that are ignored 
in this study. Although multiple studies have shown fabrics have a no¬ 
ticeable effect on the overall drag of a human body [7, 38], the walls in 
this computational work can be regarded as smooth. In all simulations, 
the models and ground of the domain are modeled as non-moving walls 
with no-slip conditions. The clothing that is part of the models is treated 
in the same manner. 


Table 1. Summary of boundary conditions in the present study. 


Wind speed 



m-s 1 

1.0, 2.5, 5.0, 7.5, 

10.0, 15.0, 20.0, 

25.0, 30.0 

Reynolds 

Re 


- 

1.281 ■ 10 5 ~ 

Number 




3.384 • 10 6 

Domain bounds 

— 


m 

r-1.44 < x < 1.44 
< -2.34 < y < 5.65 





l 0 < z < 2.32 

Turbulent 

intensity 

frnlet 



1% 


^outlet,backflow 

- 

5% 

Turbulent 

/aO 

\ 

- 

10 

viscosity ratio 

\»J 

/aO 

inlet 

) 

- 

10 


\AU 

outlet,backflow 



Outlet gauge 
pressure 

p g 


Pa 

0 


Inlet velocities range from 1 m-s' 1 to 30 m-s' 1 in the positive y- 
direction (since in this respect, the positive z-direction refers to the “up” 
orientation; refer to Figure 4 for clarification), highlighting typical wind 
speeds encountered on a day-to-day basis such as walking [39] all the 
way up to standing in a violent storm [40]. At the inlet, turbulence is 
specified using both turbulence intensity, /, and turbulent viscosity 
ratio, g t /g. Turbulence intensity is defined as the ratio of the root-mean- 
square of velocity fluctuations, u', to the mean flow velocity, U^, and 
the turbulent viscosity ratio being directly proportional to the turbulent 
Reynolds number (Re t = k 2 /ev). These values are summarized in Table 
1 . 

The boundary condition at the outlet is treated as a pressure outlet 
where a static gauge pressure is specified. In this case, turbulence is 
specified similarly as the inlet but regarded in terms of “backflow”, 



Figure 5. Side view of the full grid domain along the median plane. 


should the flow reverse direction at the boundary during iterative calcu¬ 
lations. The remaining borders of the “virtual wind tunnel” are modeled 
as symmetric to simulate zero-shear slip walls. In FLUENT, this bound¬ 
ary condition assumes a zero flux for all quantities, which imposes a 
zero normal gradient across the defined boundary and thus enforces a 
parallel flow. 

In FLUENT, the flow is initialized with a velocity field equal to the 
specified velocity for the run, e.g., a run specified at 1.0 m-s' 1 would 
have the entire field initialized with that value, and so on. Turbulence 
parameters at the boundaries are also initialized based on turbulence 
values as specified in Table 1 . The blockage ratio was determined to be 
8.7%, which would necessitate the usage of a correction factor to data; 
however, a blockage ratio of up to 10% in regards to bluff bodies has 
shown to provide reasonably similar outcomes compared to testing 
using lower blockage ratios [41] and therefore a correction factor was 
not used. 

3.2 Grid Generation 

The computational domain was discretized with an unstructured 
grid as shown in Figure 5. To reduce numerical diffusion and to more 
accurately resolve the viscous boundary layer, the surface grids on the 
models and ground were extruded using prismatic elements that are 
sized appropriately to the aspect ratio of their associated surface cell. 
These prisms are grown to 5 layers and follows recommendations put 
forth by Lanfrit [42]. Two prismatic bodies of influence (BOI) of in¬ 
creasing refinement are used to improve the resolution of the grid in 
both the wake region and the surrounding area around the models to 
sufficiently capture turbulence and flow separation. This is done to en¬ 
sure that computational processing is focused on more important re¬ 
gions in the flow regime while keeping the far field sufficiently coarse 
enough as to not dramatically hamper computational time. The overall 
grid is limited to a maximum spacing of 0.1 m and a minimum of 2 mm. 
The smaller, finer BOI is sized by the bounds -0.894 < x < 0.894, 
-0.1 < y < 2.76, 0 < z < 1.87 and the larger, coarser BOI defined by 
-0.944 < x < 0.944, -1 < y < 4.26, 0 < z < 1.92, all in meters. 

A conversion algorithm in FLUENT was used to convert the pre¬ 
liminary tetrahedral and prismatic grid into a polyhedral one. Polyhedra 
exhibit advantages over tetrahedra, namely, they approximate gradients 
better than tetrahedra due to the fact they are bounded by many neigh¬ 
bors. Additionally, polyhedra have more lax geometric criteria due to 
their insensitivity to stretching, making grid pre- and post-processing 
easier; this is well suited to the highly complex geometry of the models 
used. It has been observed that polyhedral grids provide the same level 
of accuracy as tetrahedral ones, but of significantly lower element 
count, thereby hastening simulations [43]. Furthermore, polyhedral 
grids have shown to improve convergence while having notably greater 
accuracy under unsteady simulations [44]. This is further supported by 
similar external aerodynamic studies run under FLUENT, where 
speedups between 2 to 3 times towards a converged solution have been 
observed [45]. 

3.3 Turbulence Model and Computational Approach 

The flow around the models is modeled with Reynolds-averaged 
Navier-Stokes (RANS) equations in incompressible form. Written in 
Cartesian tensor form and having flow variables of the form </) = </> + <// 
(with cp and (//being the mean and fluctuating components respectively) 
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being substituted into the instantaneous continuity and momentum 
equations of the exact Navier-Stokes equations, we obtain [46] 


^ + 4-0>“;) = ° 

dt dX; 
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The closure of the set of equations are done by means of Menter’s 
two-equation blended k-s / k-co shear stress transport (SST) model [47] 
which computes eddy viscosity with a linear stress-strain closure. Ment¬ 
er’s model combines and smoothly blends the individual strengths of 
the k-s and k-co models, where the k-s model sufficiently predicts turbu¬ 
lence in both the free stream and wake, and the k-co model more accu¬ 
rately predicts boundary layer separation in adverse pressure gradients 
near no-slip surfaces. 

The spatial convective terms in Equations (9) and (10) are discre¬ 
tized using a second-order upwind scheme (except for pressure being 
solved using PRESTO! — PREssure STaggering Option, which is simi¬ 
lar to staggered-grid schemes used on structured grids [48]) with the 
diffusive terms discretized using a weighted least squares cell-based 
construction technique. The inadequacies of the least squares approach 
are well-known; however, it provides accuracy that is comparable to 
nodal schemes [49] and due to the use of a polyhedral mesh, the least 
squares approach provides an adequate balance between computational 
expense and accuracy. The solver utilizes a pressure-based coupled 
approach to the equations, which is found to have superior convergence 
and reduced computational time than segregated algorithms [50]. The 
computational fluid dynamics code used in ANSYS FLUENT has been 
validated and used in many widespread applications, e.g., reverse flow 
in converging channels [51], moment and lift predictions on full-scale 
3-D models of airplanes [45], bluff-body drag prediction [52], heat 
transfer in Couette flows [53], and complete direct injection internal 
combustion engine simulation [54]. 

Simulating the flow around a human body was found to be naturally 
unsteady as shown by Edge et al. [55], and the ideal approach would be 
solving the flow equations in a time-dependent fashion then finding the 
computed mean quantities manually. However, simulating the flows in a 
time-accurate manner and finding the computed means would be pro¬ 
hibitively time-consuming for this study, especially since there are a 
total of 18 different runs that would need to be completed. In addition to 
this, obtaining time-averaged results from transient simulations for all 
of the required runs would require computational resources beyond the 
capabilities of the author’s reach. Thus, the most realistic and only 
feasible approach was to attempt a steady-state calculation, then deter¬ 
mine iteration-averaged quantities once a quasi-steady-state in the flow 
field was developed. This method is not without major drawbacks; no¬ 
tably, obtaining exact quantities can only be achieved by finding the 
limit of an infinite sample and it is therefore expected that the comput- 
ed-mean will contain sampling errors. Moreover, this approach would 
more than likely imply that small shedding features are not resolved and 
thus their effects on the overall flow are ignored; this suggests that there 
will be accuracy implications in terms of solutions to the problems un¬ 
der examination. 

In order to give the simulations a semblance of a fighting chance, a 
pseudo-transient time-stepping method was applied to the flow equa¬ 
tions being solved. This implicit under-relaxation serves as a predictor- 
corrector method with the objective of fast convergence and accurate 
temporal integration that achieves an approximate steady state [56]. 
Furthermore, the use of pseudo-transient under-relaxation factors 
(URFs) was found to accelerate convergence in a more robust manner 
[57] and the rate of convergence has been shown to reassuringly arrive 
at a solution earlier than without the use of these URFs [58]. Hence, by 
taking advantage of the natural structure of the problem by evolving the 


5 The astute reader should note that for all of the simulations, a laptop using 
2011-era hardware (Intel i7-2820QM processor, 16 GB RAM) was used. 
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Iterations 

Figure 6. Plot of drag coefficient versus number of iterations on a medium mesh 
with an inlet velocity of 10 m-s' 1 . 


flow equations in time, we are able to achieve a solution that reasonably 
integrates transient effects in a way that can be realized through the use 
of a steady-state calculation. In FLUENT, the pseudo time step was left 
to be calculated automatically following the recipe found in Section 
21.6.1 of the Theory Guide [46]. This is done to hasten setup and retain 
consistency across all 18 simulations. 

3.4 Convergence Discussion 

Initial calculations were run to preview the exact behavior of the it¬ 
erative calculation process; indeed, as previously demonstrated by Edge 
et al., the inherent unsteadiness of the flow resulted in difficulties ob¬ 
taining reduced residuals, with the scaled continuity residuals hovering 
below but near 1-10' 2 , and other residuals remaining below 1-10' 4 . A 
technique stemming from experience was devised to assist in reducing 
the residuals while maintaining spatial accuracy. For the first 50 itera¬ 
tions, a blending method between first and second order spatial discreti¬ 
zation schemes (with a bias for the first order scheme), lower URFs, and 
an “aggressive” pseudo-time step, were used to rapidly approach a qua¬ 
si-stable solution. The following 100 iterations switched the spatial 
discretization scheme to a higher order method along with raised URFs 
to increase accuracy. The remaining iterations were then run using a 
“conservative” pseudo-time step scaled by 0.5 to further converge and 
resolve flow details that are more apparent on a smaller timescale. This 
technique enabled the residuals to drop by nearly an order of magnitude, 
and thus it was used for all other simulations. 

As discussed in Section 3.3, taking mean quantities of a finite sam¬ 
ple will incur sampling errors of which are difficult to quantify com¬ 
pared to a mean derived from an infinite sample. However, a potential 
approach to help determine the exact nature of these errors would be 
taking an integrated surface quantity, such as drag coefficient, and see¬ 
ing how this quantity changes as the iterations progress. Such a measure 
serves as an observed global variable in that the measured quantity is 
dependent on how accurate the grid can resolve the flow. Figure 6 
shows the iteration history of drag coefficient between the two models 
taken with an inlet velocity of 10 m-s' 1 and run using a grid with medi¬ 
um refinement (refer to Section 3.5 for refinement details). The figure 
shows that the solution reaches quasi steady-state values after 150 itera¬ 
tions, with an oscillatory period occurring roughly every 140-150 itera¬ 
tions. A 100 iteration moving average calculated by FLUENT was ap¬ 
plied to the values to assist in determining if overall convergence was 
reached. Looking at the tendency of the averages, it is reasonable to 
conclude that the solution is indeed quasi-stable. With the average vary¬ 
ing very little compared to actual iteration values, and for the sake of 
minimizing computational expense in this study, this observed tendency 
incurs a reasonable assurance that the simulations are adequately con¬ 
verged. Finally, the net flux imbalance was calculated to be less than 
1 % of the smallest flux through the domain, further supporting conver¬ 
gence [59]. 
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3.5 Numerical Uncertainty and Grid Selection 

Moving forward, to preserve accuracy and to assure low computa¬ 
tional costs, the effect of grid resolution is studied by comparing results 
on drag coefficient using different grids of different refinement. The 
importance of grid convergence and examination of discretization errors 
in CFD simulations have been demonstrated across numerous organiza¬ 
tions [60, 61] and journals [62], and the necessity of quantifying these 
errors have led to well-established methods that attempt to describe 
effects of resolution to that of the extrapolated solution [63]. Although it 
has been repeated in this paper that quantifying exact errors is essential¬ 
ly a fruitless adventure, the additional time spent on understanding the 
potential discretization error nonetheless is useful in helping select a 
grid that balances the required resolution to minimal processing power. 

In this grid refinement study, each model used three different grids 
that are refined methodologically by decreasing the cell spacing within 
the BOIs. Between both models, the grid generation methods remained 
the same and follows the procedures as described in Section 3.2. Repre¬ 
senting fine (grid 1), medium (grid 2), and coarse (grid 3) respectively, 
the Normal model utilized grids containing 2.82-10 6 , 2.32-10 6 , and 
1.53• 10 6 elements. The Flat model used 2.67-10 6 , 2.16* 10 6 , and 
1.41 * 10 6 , elements. The difference in overall element counts are at¬ 
tributed to the reduced surface area of the Flat model compared to the 
Normal one. For sake of completeness and following the method as 
recommended by The American Society of Mechanical Engineers [64], 
the subsequent equations are presented to determine the discretization 
error between the described meshes. Letting h denote the representative 
grid size, 


h = 


N 



( 11 ) 


where AFJ is the volume of each cell. Let h x < h 2 < h 3 and let r desig¬ 
nate the refinement ratio between successive grids, such that r 21 = 
h 2 /h l and r 32 = h 3 /h 2 . Calculating the apparent order of accuracy, p , 
involves solving the following expression using a fixed-point iteration 
and using the first term as the initial guess: 


P = 


1 

ln(r 2 i) 


\n\e 32 /£ 2l \+ln 



( 12 ) 


determine a sufficiently averaged quantity. The results from the previ¬ 
ous equations are summarized in Table 2. 


Table 2. Grid summary and calculations of discretization error between the two 

models. 


4> = drag coefficient at 

10 m-s -1 

Lucoa, Normal 

Lucoa, Flat 

n u n 2 ,n 3 

2.82 • 10 6 , 2.32 • 10 6 , 

2.67- 10 6 , 2.16- 10 6 , 


1.53 • 10 6 

1.41 • 10 6 

Average wall clock 
time per 100 iterations 

72,38, 19 

67, 34, 17 

(minutes) 

Average memory usage 
per compute thread 

1674, 985, 648 

1546,916, 542 

(MB) 
h\,h 2 , h 3 

0.025541,0.027221, 

0.026007, 0.027912, 


0.031327 

0.032206 

01» 02’ 03 

0.9352,0.9418, 1.0015 

0.9501,0.9581, 1.0178 

r 2l> r 32 

1.06577, 1.15083 

1.07325, 1.15384 

£ 21’ £ 32 

0.00659, 0.09700 

0.00798, 0.05974 

P 

12.8180 

11.4372 


0.9299, 0.9299 

0.9437, 0.9437 

zw 

0.70%, 6.34% 

0.84%, 6.24% 

e 21 e 32 
c ext’ e ext 

0.56%, 1.27% 

0.68%, 1.53% 

gci 21 ,gci 32 

0.70%, 1.57% 

0.84%, 1.88% 


All grids were noted to converge in a monotonic manner, indicating 
the average flow field certainly benefits from grid refinement. As shown 
by the lower values of GCI 21 to that of GCI 32 , it would be most benefi¬ 
cial to run all the simulations using the fine grid. On the other hand, as 
accuracy increases, so do memory and compute requirements. The fin¬ 
est grid managed to exceed the available physical memory on the ma¬ 
chine that the simulations were run on. In addition, calculations were 
found to require nearly double the amount of time to arrive at a solution 
compared to the medium grid. The tradeoff between the reduction in 
GCI and resource consumption was too great considering the limited 
resources to begin with and the potential amount of time required for all 
18 simulations. Therefore, as shown from this grid study, the medium 
grid serves as the most practical balance between accuracy and resource 
expenditure; for all simulations, the medium grid was selected. 


where s = sign(£ 32 /e 21 ), e 32 = 0 3 - </> 2 , e 2l = <p 2 - <p x , and cp i repre¬ 
senting the quantity being investigated. Finding the extrapolated, as¬ 
ymptotic value involves calculating 

4>? x , = (r P 2l 'P l -4> 2 )/(S 2l -l) (13) 

with 4> 3 e l t being similar. Then, the following error estimates can be 
found: 


4 >\ — 


0 i 


0&-0i 


GCIo 


€a 

A - 1 


(14) 

(15) 


(16) 


4. Results and Discussion 

The immediate discrepancies between Lucoa’s model and other 
human models used in similar studies (other than the fact that she has an 
enormous chest) are the presence of horns, a baseball cap, raised arms, 
and solidly modeled hair. Each affects the wake in their respective man¬ 
ner, with these individual effects being described in the following sub¬ 
sections. It is important to recognize that although these differences 
between Lucoa and other human models are present, they nonetheless 
do not affect the direct comparison between the Normal and Flat mod¬ 
els, which was a primary objective of this work. 

The following subsection details the variances of drag and lift as air 
velocity changes and compares the outcomes to previous studies on the 
human form. The subsequent subsections delve into the differences in 
wake structure through the analysis of streamlines, velocities, TKE, and 
skin friction. While no sound comments can be made on the concise 
form of instantaneous flow structures, quasi steady-state general fea¬ 
tures present in the flow are described. 


with e 2 Q describing the approximate relative error, e 2 ^ xt describing the 
extrapolated relative error, and GCI 21 being the grid convergence index 
(GCI). e 32 , e 22 xV and GCI 32 are calculated similarly. 

The GCI is an indicator with a 95% confidence interval of how far 
the finer of the two compared grids is to the asymptotic value (p ext and 
predicts how further refinement affects the solution [65]. The evaluated 
quantities are taken after 400 iterations, which in Figure 6 and discussed 
previously, has shown to provide a reasonable amount of iterations to 


4.1 Velocity-Based Trends 

The convergence method presented in Section 3.4 was used for all 
velocities listed in Table 1 between the two models. The drag and lift 
were calculated using the same 100 iteration moving average in 
FLUENT and the values at the end of 400 iterations were used. Figures 
7-9 plot the results between the two models. In all of the figures, a cubic 
spline was used to interpolate the data. The figures labeled with “(a)” 
correspond to drag-related values while figures labeled with “(b)” cor¬ 
respond to lift-related values. 
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Re xlO 6 Re xlO 6 

(a) Drag coefficient vs. air velocity. (b) Lift coefficient vs. air velocity. 

Figure 7. Plots comparing coefficient quantities between the two models. 



Re 


x 10 


Re 


(a) Drag force vs. air velocity. (b) Lift force vs. air velocity. 

Figure 8. Plots comparing the forces experienced on the two models at various air velocities. 


xlO 6 



Re xlO 6 

(a) Relative drag difference of the Flat model to the Normal model. 


Re 

(b) Relative lift difference of the Flat model to the Normal model. 


x 10 


Figure 9. Plots indicating the relative force differences of the Flat model compared to the Normal model. 


Figure 7 focuses on the coefficients, and it is readily apparent that 
at all velocities, the Flat model incurs higher coefficients than the Nor¬ 
mal model. In the region near 5 m-s' 1 , the drag coefficients between the 
two models were closest to one another. Judging from the errors de¬ 
scribed in Section 3.5, it is safe to presume that these differences are 
within the error bounds and therefore the difference between the two 
models near this velocity is negligible. In 7(a), it is interesting to ob¬ 
serve that the drag coefficients decrease following a power relation 
corresponding to velocity. This relation most likely arises from the pe¬ 
culiarity of the solidly modeled hair; downwash due to the hair diverts 
air behind the torso then directs the air downwards, which causes less 
air to flow directly behind the body. Presumably, this downwash closes 
off the size of the primary recirculation region (PRR) behind the body. 


As the air velocity increases, the effect of the downwash becomes more 
prominent and therefore a general decrease in the drag coefficient is 
observed. The large breasts on the Normal model may also play a role 
in enhancing this downwash effect, as the air is gradually diverted 
around the bust and is able to maintain momentum before being redi¬ 
rected by the hair. The stronger resulting downward jet of air thus closes 
off the recirculation region more effectively. The degree of which this 
downwash affects drag is difficult to determine from this study, but it is 
reasonable to deduce that due to the drag reduction of the Normal model 
compared to the Flat model, the effects are markedly noticeable. The 
geometrical presence of the breasts may also have an effect in reducing 
drag as they smoothly redirect a portion of the air around the torso, 
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whereas the Flat model simply forces majority of the air to stagnate at 
the chest, which incurs higher drag. 

The difference in lift is notably more dramatic than drag. In Figure 
7(b), the lift coefficient for both models gradually increases from 1 to 
10 m-s' 1 then levels off as the velocity changes. This increase can also 
be contributed to the downwash from the hair, since as the velocity 
increases, there is a higher flow rate and thus more fluid mass pushing 
upwards against the models. The brim of the baseball cap on both mod¬ 
els may additionally contribute to lift as it lies orthogonal to the mean 
flow direction and causes a buildup in pressure in front of the face. The 
Normal model provides a lower lift coefficient possibly due to the 
breasts providing downforce; since the Normal model has a significant 
and primarily upward sloping surface in front of the torso that the Flat 
model lacks, the downward force from the air flowing atop the bust 
most likely provides this reduction in lift. Granted, due to the geometry 
being modeled as completely solid, this effect does not account for 
breast deflection caused by these forces. Nonetheless, the presence of 
Lucoa’s large breasts evidently provides a reduction in lift. 

Figure 8 provides a comparison of the actual forces experienced on 
the models. As seen in both 8(a) and 8(b), the forces increase following 
a roughly squared relation, which is what would be expected looking at 
the velocity relationship in Equations (1) and (5). As elaborated in the 
previous paragraphs, the difference in drag is less dramatic than that of 
lift, and this is reflected by the more difficult to discern plots in 8(a) 
than in 8(b). This observation is further expanded upon in Figure 9, 
which compares the relative force difference the Flat model encounters 
compared to the Normal model. Looking at 9(a), the Flat model experi¬ 
ences a nearly 4% increase in drag at 1 m-s' 1 , however, this value de¬ 
creases immediately and falls to below 1% at 5 m-s' 1 . From 5 to 10 
m-s' 1 , the difference between the two widens again, levels off from 10 
to 20 m-s' 1 , then impetuously increases at 25 m-s' 1 and remains marked¬ 
ly greater at the highest tested velocity of 30 m-s' 1 . The reason as to 
why the drag difference reaches its minima at 5 m-s' 1 cannot be directly 
inferred from this work. Throughout the tested velocity range, the Flat 
model incurs an average drag increase of 2%, which, in the context of 
day-to-day life at low air velocities, is practically insignificant. Howev¬ 
er, in regards to performance, the large breasts on the Normal model can 
prove to be advantageous, assuming the overall position remains similar 
to the erect postures of the models. The degree of which this advantage 
can be quantified depends heavily on the type and duration of competi¬ 
tion to be examined, which is not discussed further in detail. The rela¬ 
tive lift force difference averages 24%, with the values ranging from 
32% at 1 m-s' 1 , gradually decreasing to 21% near 15 m-s' 1 , and then 
slightly increasing to 23% at the highest tested velocity. From this re¬ 
sult, the presumed effect of Lucoa’s breasts reducing lift is more pro¬ 
nounced at lower air velocities, while still remaining effective as the 
velocity increases. 


Table 3. Comparisons of drag coefficient with other researchers’ results. 


Researchers 

Year 

Reference 

Posture 

C D 

This study 

- 

- 

Standing 
(Normal model) 

0.89 ~ 1.05 




Standing 
(Flat model) 

0.91 ~ 1.09 

Hill 

1927 

[i] 

Standing 

0.98 

Schmitt 

1954 

[V] 

Standing 

(clothed) 

Approx. 1.36 




Standing 

(nude) 

Approx. 1.20 

Shanebrook and 

1976 

[66] 

Running (cylin¬ 

1.2 

Jaszczak 



drical model) 


Penwarden et al. 

1977 

[67] 

Standing 

(clothed) 

1.18 

Davies 

1980 

[2] 

Running 

0.87 

Brownlie et al. 

1987 

[68] 

Standing 

0.96 ~ 0.98 

Gomez et al. 

2013 

[69] 

Running 

1.20 

Inoue et al. 

2016 

[70] 

Running (without 

1.17 


ground effects) 


Table 3 summarizes the encountered drag that both models experi¬ 
enced throughout the tested velocity range, together with other re¬ 
searchers’ experiments. As a result, the present study reveals that Lu¬ 
coa’s form tends to correspond favorably to experiments done with 
standing models. It is seen that Lucoa’s results overlap with those of 
Hill [1] and Brownlie et al. [68]. Schmitt’s [7] results were an atypical 
case as the drag coefficients in that study were based on the entire sur¬ 
face area of the body divided by the product of volume and subject 
height, instead of the customary frontal projected area used elsewhere. 
This resulted in coefficients ranging from 10 to 13, necessitating a re¬ 
calculation in order to more directly compare to the other results in 
Table 3. It was demonstrated through Schmitt’s results that the use of 
nude subjects (akin to having smooth walls as mentioned in Section 3.1) 
result in lower drag. Although the values for nude subjects were still 
greater than those found in the current study, this finding provides addi¬ 
tional support as to why Lucoa’s overall form tends to encounter lower 
coefficients in contrast to other studies. With the objective of obtaining 
C D values for a wide range of people in everyday situations, work done 
by Penwarden et al. [67] confirms that clothing certainly increases drag. 
Furthermore, the authors in that study comment that extrapolating their 
data to simulate bare-skinned models would result in values more akin 
to that of Hill’s data, thereby supporting the notion that Lucoa’s results 
are indeed realistic. 

Comparisons to results based around running are also provided, as 
this form of competition is a natural progression from simply standing. 
Understandably, a runner will experience higher drag due to their ever- 
changing position and the strict anti-symmetry in orientation. Shane- 
brook & Jaszczak [66] investigated drag through the use of a general¬ 
ized cylindrical human model, analogous to the assumptions described 
in Section 1.1. As such, their results parallel those of other investiga¬ 
tions done on runners and show that a cylindrical model can provide 
convincing data. Davies’ [2] work discovered that the effects of Reyn¬ 
olds number remained constant below velocities of approximately 18 
m-s' 1 . Conversely, this effect was not seen in this study, which most 
likely harkens back to the drawback of the modeled hair. Gomez et al. 
[69] developed a parametric model of drag that considered both veloci¬ 
ty, w, and acceleration, w 2 , based on the record-setting performance of 
Usain Bolt. They revealed that 92% of energy used in running is ab¬ 
sorbed by drag, which bolsters both the significance of the drag reduc¬ 
tion and the implicational anthropomorphic advantage that Lucoa has. 
Inoue et al. [70] examined the effects of a moving belt system to simu¬ 
late ground effects of a runner, in addition to providing data on how leg 
position (e.g., forward or behind) changes drag. From their study, it was 
observed that having the legs placed either forward and behind veritably 
increases drag due to the posture. Consequently, it is expected that a 
prospective study incorporating the effects of running based on the pre¬ 
sent results (as described in herein and Table 3) would feature compara¬ 
ble C D values. While the current stationary results may not be directly 
applicable to the moving case, they nevertheless offer a suitable starting 
point on which to base further research. 

4.2 Salient Flow Analysis 

Discussion pertaining to the variances in flow structures is exam¬ 
ined under an air velocity of 10 m-s' 1 . In addition to being previously 
examined in Sections 3.4 and 3.5, this velocity is readily achievable 
under human locomotion [69], making the selection straightforward. 
While the results are derived from a quasi-steady-state approach, gen¬ 
eral flow structures can be reasonably interpreted. 

4.2.1 Streamlines and Flow Velocities 

As reflected in Figures 10-13, the fluid dynamics of the wake are 
shown with the use of streamlines colored according to velocity magni¬ 
tude. The figures labeled as “( a )” and “(b)” correspond to the Normal 
and Flat models respectively. The overall wake structures revealed in 
Figure 10 fall into two distinct regions based on height, with the PRR 
found above the hips and extending to the top of the head, and a second, 
more complex regime behind the legs. The works described in Section 

1.1 correctly predicted these structures. Specifically, the PRR for both 
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(a) Wake structures associated with the Normal model. (b) Wake structures associated with the Flat model. 

Figure 10. Streamlines emanating at x = 0 and spanning a line from 0 < z < 2.3 meters showing primary flow structures. Pressure experienced on the models is also 

shown. 




(a) Set of vertically aligned streamlines on the Normal model. (b) Set of vertically aligned streamlines on the Flat model. 

Figure 11. Detailed view of flow structures around the chest of the models using the same streamlines from Figure 10. 



(b) Set of horizontal streamlines on the Flat model at the same z-location. 
emanating at z = 1.18 meters and spanning a line from -0.5 < x < 0.5 meters. 


(a) Set of horizontal streamlines on the Normal model flowing atop the bust. 
Figure 12. Detailed view of flow structures around the chest using streamlines 



(a) Set of horizontal streamlines on the Normal model flowing beneath the bust. (b) Set of horizontal streamlines on the Flat model at the same z-location. 


Figure 13. Detailed view of flow structures around the chest using streamlines emanating at z = 1.14 meters and spanning a line from -0.5 < x < 0.5 meters. 
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models was found to extend roughly 1.4 meters downstream, or approx¬ 
imately 4 times the width of the torso, which agrees strongly with Oka- 
moto’s and Sunabashiri’s findings. The type of vortex shedding in this 
recirculation region cannot be determined, but it is expected to nonethe¬ 
less be asymmetric, with the large-scale structures being advected 
downstream based on the findings of Edge et al. The wake behavior 
behind the legs were also discovered to obey the findings summarized 
by Zhou and Alam, with the region from the ankles to the knees having 
wake structures associated with each leg, and the region from the knees 
to the hips behaving as a unified bluff-body. The figure also shows a 
distinct vertical “jet sheet” stemming from the gap between the legs, 
with an average magnitude approximately 35% greater than the free 
stream velocity. A vortex pair originating from the top of the legs is also 
evident. This vortex pair can be seen being caused by the downwash of 
the hair; the momentum of the downward flow interacts with the air 
flowing past the top of the thighs, which results in these prominent 
structures. Additional momentum provided by the jet of air between the 
legs further enhances the strength of these vortices. Furthermore, the 
effects of the downwash are evident in the way the vortex pair assumes 
a downward angle as the flow advects; this is indicative of the signifi¬ 
cant amount of air being captured by the geometry of the hair. 

From what can be drawn from comparing Figure 10(a) and 10(b), 
the Flat model has a weaker, less organized vortex stemming from the 
legs. This observation may further explain why the Normal model has 
lower lift than the Flat model. Aider et al. [71] described the effect of 
vortex pairs affecting lift and drag, and shown that inflow caused by 
streamwise counter-rotating vortices result in a net downward force. 
Thus, the weaker vortex pair on the Flat model contributes less towards 
reducing lift than the Normal model. Following this, it is reasonable to 
suggest that the behavior of the drag curve in Figure 7(a) can be at¬ 
tributed to the formation of these vortices; the energy being redirected 
into the formation of these structures diverts a portion of the mean flow 
away from the PRR, leading to drag reduction as the velocity increases. 

Figure 11 provides a close-up look at the same streamlines in Figure 
10 around the torso of both models. As can be seen in 11(a) and correct¬ 
ly determined in Section 4.1, the breasts provide a gradual interface for 
the air to move around the torso compared to the relatively abrupt ob¬ 
struction the Flat model imposes as reflected in 11(b). This is further 
exemplified in Figure 12, which shows a set of horizontal streamlines 
positioned at z=1.18m. In 12(a), the air flows smoothly atop the breasts 
at a velocity approximately 40% to 60% of the free stream velocity, 
whereas the Flat model in 12(b) simply forces a relatively larger portion 
of the incoming air to stall at the torso. It is interesting to note that a 
small recirculation region develops directly atop the breasts, which most 
likely is dependent on the angle between them and the chest. From what 
can be seen in 12(a), this small recirculation region prevents a small 
portion of air from surmounting the shoulders and instead redirects the 
air downwards and off to the sides of the torso. Figure 13 provides an 
upwards facing view of the models along with horizontal streamlines 
emanating from z=1.14m. In 13(a), air flowing beneath the bust gains a 
velocity that is approximately 20% higher than the mean free stream 
velocity before being diverted perpendicular to the body and interacting 
with the incoming air stream and the rest of the torso. By the action of 
the breasts, the air is diverted in such a way that it is able to maintain 
momentum while it moves around the rest of the chest. This is in con¬ 
trast to 13(b), where a portion of the air is forced to stall in front of the 
chest before being redirected around the torso, much like as was seen in 
other figures involving streamlines with the Flat model. 

Figure 14 compares the velocity magnitude of the wake behind the 
two models at five different z-locations. At these locations, data is sam¬ 
pled along a line centered at the midline (x=0) and spans downstream 
from -1 < y < 2 meters. The z-locations are taken at heights of z = 
0.1m, z = 0.5m, z = 0.9m, z = 1.18m, and z = 1.55m, all of which corre¬ 
spond to the feet, legs, hips, chest, and head, respectively. From this 
comparison, the differences in wake velocities are slight, with the only 
significant difference occurring at the height of the legs. At this height, 
it can be seen that the Flat model has a wake velocity that is consistently 
higher than that of the Normal model. This observation can be attributed 
by the lower energy present in the vortices generated by the Flat model; 



Figure 14. Comparison of wake velocity magnitudes at 5 z-locations centered 
along x = 0 and ranging from -1 < y <2 meters. Lucoa positioned and to scale 

with x-axis. 


Feet z = 0.1m, Legs z = 0.5m, Hips z = 0.9m, Chest z = 1.18m, Head z = 1.55m. 

these weaker formations are affected by the mean flow to a greater de¬ 
gree and thus the velocity magnitude in this region is higher than if the 
vortices were to have greater strength, as seen with the Normal model. 

4.2.2 Turbulence Kinetic Energy 

Figure 15 shows a slice of the domain through the median plane of 
the models, indicating the TKE present in the flow. At first glance, con¬ 
tours for both models present the same general structures as described 
earlier, such as the PRR behind the torso and the smaller structures as¬ 
sociated with the legs. Due to the automatic scaling of the colors, the 
differences present in the PRR of the Flat model compared to the Nor¬ 
mal model are slightly exaggerated. However, the differences are still 
distinguishable in that the PRR for the Flat model has a measurably 
higher level of TKE than that of the Normal model. Indeed, 15(b) indi¬ 
cates that the two apparent bands of significant TKE associated with the 
PRR are generally more intense, even after factoring the differences 
between the color scales. Common to both models is a region behind 
the hips where the relatively highest TKE was measured. This lower 
band within the PRR was found to be more compact for the Normal 
model, whereas the Flat model had a marginally less intense and more 
widespread band. This effect may be attributed to the way the down- 
wash from the hair interacts with the PRR. 

Directly below the PRR is the TKE present in the vortex pairs orig¬ 
inating from the legs. In 15(a), the Normal model has a notably higher 
level of TKE associated with these vortices, coupled with the observa¬ 
tion that this energy is maintained as the flow advects downstream. In 
contrast, 15(b) shows that the TKE in the vortices is lower and that the 
energy is dissipated sooner. This remark further supports the interpreta¬ 
tions hitherto; that the formation of stronger vortices the Normal model 
generates plays a role in reducing both drag and lift by diverting energy 
away from the PRR. 

Smaller regions with notable TKE present between both models are 
those associated directly behind the legs, at the ankles, and directly 
above the brim of the baseball cap. Between the two, the TKE at the 
ankles remains unperturbed by the differences in the wake structures 
above it. This region is related to the individual wakes associated with 
each foot, behaving much like the two cylinder arrangement as de¬ 
scribed earlier. A downward “jet” of TKE originating from the thighs 
and curling parallel to the mean flow direction at the height of the knees 
shows some variability between the two models. This region is linked to 
the “jet” of air that is formed between the legs interacting with the 
downwash from the hair. With the Normal model in 15(a), this “jet” of 
turbulence remains both vertical and closer to the legs. In 15(b), this 
region takes on a slightly more horizontal transition and extends farther 
into the wake. A small recirculation region, much like that associated 
with the breasts on the Normal model, is found above the brim of the 
baseball cap, which is due to the blunt transition the shape of the head 
imposes to the incoming flow. 
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(a) TKE associated with the Normal model. (b) TKE associated with the Flat model. 

Figure 15. TKE present in the flow behind the models along the median plane. 



(a) TKE present behind the breasts of the Normal model. (b) TKE present around the torso of the Flat model. 

Figure 16. TKE present around the torso of both models along the coronal plane. 
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Figure 17. Comparison of TKE at 5 z-locations centered along x = 0 and ranging 
from -1 < y <2 meters. Lucoa positioned and to scale with x-axis. 

Feet z = 0.1m, Legs z = 0.5m, Hips z = 0.9m, Chest z = 1.18m, Head z = 1.55m. 

Figure 16 shows the TKE around the torso of both models along the 
coronal plane. As obvious as the difference between both models, the 
presence of a turbulent region behind the breasts and below the axilla 
comes as no surprise. It is apparent that this region is influenced by the 
small recirculation region that forms above the breasts spilling air 
downwards and behind the torso, which can be seen in Figure 12(a). 
Furthermore, air from beneath the breasts also contributes to the size of 
this region, as in Figure 13(a) it can be understood that this layer of 
accelerated air provides additional shear against the flow from above 
and increases the turbulence behind the breasts. Since the Flat model 
lacks such features, Figure 16(b) presents an understandably unremark¬ 
able result, with only a small sliver of turbulence along the side of the 
torso. Common to both models however, is turbulence being generated 
by the hair itself, stemming from the hair being left modeled as solid. 
The true extent of how much turbulence is generated from the solid hair 
and its effects on the PRR is unfortunately intractable from this current 


study due to its highly complex geometry. Though, it is expected that 
simulations dealing with physically accurate hair may undoubtedly 
increase turbulence, and hence, drag, based on casual observations done 
on how hairstyles can affect aerodynamic performance [72]. 

The same approach used to obtain information as seen in Figure 14 
was used in Figure 17, with TKE being measured instead of velocity. As 
indicated from this figure, it is seen that the vortices originating from 
the legs indubitably diverts energy away from the PRR and supports the 
findings stemming from Figure 15. This is reflected with stronger over¬ 
all turbulence present at z-locations corresponding with the hips, chest, 
and head for the Flat model, along with markedly lower TKE at the 
level of the legs in contrast to the Normal model. Slight differences are 
noted, with the Flat model having a higher peak TKE at the legs imme¬ 
diately behind the body than the Normal model. Then, from 1 to 2 me¬ 
ters downstream, the Normal model generates a greater amount of TKE. 
This observation may ultimately be attributed to the difference in shape 
the leg vortices assume between the two models. 

4.2.3 Skin Friction Coefficient 

An evaluation of the effects of the macroscopic geometry between 
the two models on skin friction is shown with Figure 18. This figure 
represents a scatter plot of C f measured at every vertex on the models, 
with regions of relatively high C f labeled according to which location 
on the model they are found on. Two directions are evaluated, with the 
x-axis on 18(a) representing the streamwise position, and the x-axis on 
18(b) representing the z-location. A silhouette of the models present in 
the background of the plots are positioned and scaled to their respective 
x-axis to further provide context with the labeling. Additionally, the data 
was averaged and plotted along with the scatter plot to reveal discrep¬ 
ancies that may be more difficult to discern. 

Expectedly, both models feature the same general scatter plots, with 
nearly negligible differences. However, common to both models and the 
plots are the presence of numerous spikes associated with the hair. In¬ 
asmuch as they are an artifact of the hair being modeled solidly, these 
spikes are aptly denoted as “hair anomalisms” within the plots. As such, 
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Figure 18. Scatter plots of skin friction comparing both models. Sources of regionally significant C f are labeled. An additional plot denoting a smoothed average of the 
scatter data is present to provide clarity in comparison. Silhouettes of the models are positioned and scaled according to their respective axis in the plots. 


the hair is undeniably a major contributor to skin friction on the models. 
This is especially apparent in 18(a), with the source of skin friction 
being contributed solely by the hair from the rear of the model. Another 
exemplification of the hair’s effect onC f can be seen in 18(b), with a 
grouping of peaks obtaining a value of 0.06 corresponding to the hair 
tips and bangs around the face. Additional sources of significant C f can 
be easily observed in the same figure, such as the ankles, around the 
thighs, the fingertips, and especially Lucoa’s horns. These dramatic 
changes to skin friction indicate where on the models flow separation is 
more likely to occur, which in turn, relates to where turbulence can 
potentially be generated. 

The breasts, compared to the relatively noisy scatter generated by 
the rest of the body, provides a smooth transition in regards to skin fric¬ 
tion over its surface. A small spike found at the downstream location of 
y = -0.77m in 18(a) corresponds to the small recirculation region above 
the breasts seen in Figure 15(a). In a sense, it can deduced that the 
breasts act to reduce C f , which can be seen as a small dip in the 
smoothed data centered around z = 1.2m in 18(b). However, even with 
such a diverse of a plot, the average C f across the whole body can be 
seen as roughly 0.01 from both figures, which indicates how little skin 
friction actually affects overall drag with bluff bodies. 

5. Summary and Recommendations 

This paper has offered a unique compendium of data providing in¬ 
sight into the effects of a specific physiological feature on the aerody¬ 
namic performance of a human. As such, the results have indicated that 
large breasts can be notably aerodynamic through the reduction of drag 
and lift. The Flat model incurred a 4% maximum drag increase com¬ 
pared to the Normal model, with an average of approximately 2% span¬ 
ning velocities from 1 to 30 m-s' 1 . The Flat model also experienced 
more lift, with a maximum difference being 32% and averaging 22%. 
As illustrated, the mechanism behind the drag and lift behaviors ob¬ 
served between both models was elucidated through the analysis of 
streamlines around the body and the structures associated with TKE; the 
Normal model provides advantageously lower drag and lift by the gen¬ 
eration of stronger vortices from the legs, which in turn originates from 
the action of the breasts redirecting the flow around the torso. From 
what has been presented in this preliminary work, it is safe to conclude 
that the phrase “Flat is Justice” is deficient aerodynamically. 

A major shortcoming intrinsic to this study was the decision to 
leave the hair present in the models as immovable and solid. As noted 
earlier, a significant portion of the air flowing around the body was 
captured by the geometry of the hair, and this affected the wake struc¬ 
tures behind the model. This effect therefore blunts the overall applica¬ 
bility of the results found in this study to actual human models. It is 
important to note, however, that large breasts do give a consistently 
notable aerodynamic advantage, as reflected in the overall lower forces 


experienced by the Normal model. Additionally, even though the wake 
structures generated by the hair resulted in a departure from the ex¬ 
pected drag and lift behaviors, comparisons done with other experi¬ 
menters’ results show that the outcomes from this study were strongly 
promising. 

For future studies, several recommendations are provided. In the 
near-term, a re-evaluation of the current work without the hindrance of 
modeled hair should be done. Work done using a time-dependent com¬ 
putational approach should also be completed to further gauge the ef¬ 
fects and inaccuracies of using the pseudo-transient method in relation 
to drag and lift. A keynote proposition would be the use of a wind tun¬ 
nel experiment to acquire validation data. Ideally, this experiment 
would use live subjects of varying breast sizes in order to provide addi¬ 
tional data in which to compare to the data herein. Additional research 
to evaluate the degree of which the information provided from this cur¬ 
rent work applies to real world situations should also be completed. 
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Appendix A: Supplementary Material 

Additional contour plots in higher resolution and quality from the other 
velocities not examined in Section 4.2 can be found by following this 
link: https : //imgur. com/a/bz31B . 

A video further visualizing the flow structures in 3-D can be seen here: 
https://youtu.be/1414xHh6twO . 

The entire CFD study (project files and generated data) can be found by 
following: 

https://mega.nz/#!lz5hzA5Q!JpG_tOIbflLuI24sClpkV- 

C7nKW10xqZSMzBeHlMt_A . 
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Figure 19. Velocity contours and vectors of the Normal model along the median 

plane. 



Figure 20. Velocity contours and vectors of the Flat model along the median 

plane. 



Figure 21. Pressure coefficient of the flow regime around the Normal model 
along the median plane. 



Figure 22. Pressure coefficient of the flow regime around the Flat model along 
the median plane. 
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